Sensors 2014, 14, 12687-12700; doi:10.3390/sl40712687 



OPEN ACCESS 



sensors 

ISSN 1424-8220 

www.mdpi.com/journal/sensors 

Article 

Sensor Distribution Design of Travel Time Tomography 
in Explosion 

Yali Guo *, Yan Han, Liming Wang and Linmao Liu 

National Key Laboratory of Electronic Test Technology, North University of China, Taiyuan 030051, 
China; E-Mails: hanyan@nuc.edu.cn (Y.H.); wlm@nuc.edu.cn (L.W.); liulinmao@nuc.edu.cn (L.L.) 

* Author to whom correspondence should be addressed; E-Mail: guoyali@nuc.edu.cn; 
Tel.: +86-138-0343-2936; Fax: +86-0351-355-7396. 

Received: 7 May 2014; in revised form: 5 July 2014 / Accepted: 8 July 2014 / 
Published: 15 July 2014 

Abstract: Optimal sensor distribution in explosion testing is important in saving test costs 
and improving experiment efficiency. Aiming at travel time tomography in an explosion, 
an optimizing method in sensor distribution is proposed to improve the inversion stability. 
The influence factors of inversion stability are analyzed and the evaluating function on 
optimizing sensor distribution is proposed. This paper presents a sub-region and 
multi-scale cell partition method, according to the characteristics of a shock wave in an 
explosion. An adaptive escaping particle swarm optimization algorithm is employed to 
achieve the optimal sensor distribution. The experimental results demonstrate that optimal 
sensor distribution has improved both indexes and inversion stability. 

Keywords: sensor distribution design; travel time tomography; explosion; sub-region and 
multi-scale cells; adaptive escaping particle swarm optimization 



1. Introduction 

Explosion technology is applied increasingly [1,2]. As explosion experiments are complex and the 
costs are high, optimizing sensor distribution is necessary to save test costs and improve experiment 
performance. The number of sensors to be used and their locations should be optimized [3-5]. 

When a shock wave propagates in air, the shock wave velocity is very fast. Supposing that a shock 
wave propagates by direct rays without propagation by grid boundary, in the course of shock wave 
propagation, the relationship of travel time is: 
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t = j R — -dr = j R s -dr 



a) 



The wave front normal is defined as rays, i.e., r. Each sensor corresponds to a ray. v is velocity and s is 
slowness. The travel time t is the integral of slowness along with rays. 

The test region is divided into some regular or irregular cells as shown in Figure I, the 
discretization expression of Equation (1) is: 



^=T d tJ s J 0' = l,2,3,-,M;y = l,2,-,JV) 



.7=1 



(2) 



where U is the travel time that a shock wave travels from the explosive to the ith sensor, dy is the 
length of the ith ray in the jth cell. Sj is the slowness in jth cell. M is the sensor number and N is the 
cell number. Equation (2) can be written in a matrix form as: 

DS = T (3) 

where T = (t 1 ,t 2 ,---,t M )' is the m-dimension column vector of travel time; S = (s x ,s 2 ,---,s N )' is an 

unknown N dimensional column vector; D is the distance matrix of M x N and its element is dy [6,7]. 
The elements of Tcan be obtained by the tested data from sensors. The matrix D can be calculated by the 
position of the explosive and sensors. The unknown slowness S can be obtained by solving Equation (3). 
Based on the above-mentioned principle, the velocity distribution of shock wave can be obtained. 

Figure 1. Illustration of travel time tomography. 




The explosive # Sensors 

For the travel time tomography, the driving source is single and the sensors are few. The 
tomography rays are distributed sparsely. The Equation (3) is ill-posed and underdetermined [8]. For 
this tomography modality, optimizing sensor distribution is necessary for reducing costs and increasing 
acquired information. Optimizing sensor distribution is to solve ill-posed tomography problems from 
the perspective of mathematics and then the inversion stability can be ensured. 

The ill-posed degree of Equation (3) is related to the structure of matrix D. Improving the ill-posed 
degree of Equation (3) is to improve the structure of matrix D. The structure of matrix D rests with a 
parameterized model of the test region and sensor distribution [9]. 
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2. Theory of Optimal Sensor Distribution and Indexes 

2.1. Effect of Eigenvalue and Rank on Inversion 

I |2 

In Equation (3), with a given data vector t 0 , the model vector s 0 GS and then f 0 -Ds 0 \ is 
minimized. This is accomplished by pre-multiplying Equation (3) by D T and taking a matrix inverse: 

S = (D T D) l D T T (4) 

As the NxN square matrix L = D T D is often near-singular, it results in instability in the solution. 
That is, some of its eigenvectors \e t :i = l 9 -~ 9 N} have extremely small eigenvalues j/t. :/ = l,---,iV}. 

Measurement errors in data space T propagate into the solution S in parallel to each eigenvector with 
an amplification 1 / 2. . Hence, when small eigenvalues exist, the solution becomes unstable and the 

inverse problem is ill-posed [10]. When an eigenvalue is zero, the corresponding eigenvector in the 
data space can not be mapped into the model space [11]. The greater the rank is and the larger the 
eigenvalues are, the more stable the inversion problem is and the more independent pieces of 
information may be gained from the data. Hence, optimal sensor distribution can be obtained by 
maximizing the magnitude of eigenvalues of L and the rank of D. The evaluating function is: 

E x = l - — + N- rank(D) (5 \ 

trace{D T D) w 

N 

where X x is the maximal eigenvalues of D T D ; trace{D T D) = ^A t ; N is the cells number, rank(D) 

i=\ 

is the rank of matrix D . 

The evaluating function E\ has two components, one represents the relative distribution of 
eigenvalues and the other indicates the relative size of null space. If the value of E\ is minimal, the 
inversion results are optimal and stable. 

2.2. Effect of Condition Number on Inversion 

The Equation (3) is ill-posed if small perturbations of matrix DorT can result in a large change in 
solutions [12]. 

Assuming that the observed data T has a minor perturbation of ST , the perturbation of solution is 
SS . Equation (3) becomes: 

D(S + SS) = T + ST (6) 

Then, 

SS = D l ST (7) 

1 |, H 1 \\DS\\ \\T\\ 

According to the property of subordinate norm, \\SS\\ = \\D ST\\ < LD <571 and LSI > \ — rp = . It 



D D 



tell ||£~1||£r|| \\ss\\ \\st\\ 

concludes that < " 1 . That is i— I < cond(D) - 



S T / D S T 
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Then, 



cond(D) = \\DiD~ 



(8) 



where cond(D) is the condition number of matrix D. When the observed data T has a minor 

perturbation, the relative error of solution is determined by the condition number. 

Supposing that T is accurate and the matrix D has a minor perturbation of SD , the corresponding 
perturbation of solution is <SS . Supposing that S c = S + SS is the solution of perturbation equation, 

(D + 5D)-S C =(D + 5D)-(S + SS) = T ( 9 ) 
Similarly the following condition can be obtained: 

toll toll 

tt-^ — ^ < condiD)^. — ^- (]()\ 
\\S + SS\\ \\D\\ {LV) 

When the matrix D has a minor perturbation, the error of solution is also determined by the 
condition number. 

Therefore, the condition number is the second judging index of matrix D's quality. The smaller 
condition number means the more well-posed equation. 

The evaluating function expressed by the condition number can be written as: 

E 2 = cond(D) (11) 

Optimal sensor distribution can reduce the condition number and the more stable inversion 
solutions can be obtained. 

2.3. Effect of Ray Coverage on Inversion 

When designing sensor location, rays coverage should be enlarged as much as possible, i.e., rays 
density and orthogonality should be maximized. The ray density represents the number of rays passing 
through each cell. The cells not being hit by any ray are the main factors giving rise to the null space. 
They make some column vectors of matrix D be zero (dj = 0) so that the equation can not be resolved. 
Therefore, tomography with sparse rays must ensure that any column vector is non-zero. Increasing the 
ray density can avoid zero vectors. 

Ray orthogonality is measured by maximal sinusoidal quantity of angle between rays [13]. The 
small orthogonality makes some rows in D linearly dependent. 

The greater the ray density is, the better the orthogonality is and the smaller inversion error can be 
achieved. The evaluating function expressed by rays density and orthogonality can be written as: 

^-%tpj+%i.Oi-lP + k,0 ( i2) 

where p. is the ray density in the jth cell. O . is the ray orthogonality in the jth cell. The value of k x and 
k 2 are determined by the values of p. and 0. . Sensor distribution can be optimized by 
maximizing E 3 . 
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2.4. Evaluation Method 



As sensor distribution can be optimized by minimizing E x and E 2 and maximizing E 3 , an 
evaluating function is defined as: 

CO, 



E = co x E x + a> 2 E 2 +- 



(13) 



where ^ , &> 2 , and &> 3 are determined by the values of E\ 9 E2, and E3. Sensor distribution can be 

optimized in terms of minimizing E. 

There are many factors on the ill-posed and underdetermined equations except for the above analysis. 
However, the value of E can be used as a judging index with regard to optimizing sensor distribution 
and the experiment. The detailed process is as follows and a flow chart is illustrated in Figure 2. 



Figure 2. Illustration of evaluation process. 
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(1) Dividing cells according to model character. 

(2) Giving a distribution model randomly according to the number of sensors and calculating 
matrix D and the rank of D. 
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(3) If all column vectors of matrix D are non-zero and D is full rank, make this distribution 
model as an initial model; otherwise, go to Step (2). 

(4) When the number of initial model is equal to the given number K , the optimization algorithm 
is employed to obtain the optimum value of E and sensor distribution. 

3. Optimizing Sensor Distribution 

3.1. Sub-Region and Multi-Scale Cells Partition 

The test region is divided into cells and each cell has the same velocity. The dividing pattern of 
cells affects matrix D. We propose a sub-region and multi-scale cell partition method. Cells are divided 
in terms of the solution distribution, i.e., the smaller size and the higher density of cells correspond to 
the bigger changing gradient of solutions, and the bigger size and the lower density correspond to the 
smaller changing of solutions. 

When the explosive is exploding in air, the shock wave overpressure attenuates quickly in the near 
region to the explosive and the attenuation becomes slow with the increase in distance. According to 
the shock wave characteristic, the smaller cells are adopted in the near region to the explosive while 
the bigger cells are adopted in the far region. In order to avoid a row vector to be linearity dependent 
on matrix D aroused by symmetry, cells are divided into different size in the symmetrical region with 
consideration of the resolution of inversion. 

3.2. Optimizing Sensor Distribution Based on Adaptive Escaping Particle Swarm Optimization Algorithm 
3.2.1. Particle Swarm Optimization (PSO) Algorithm and Modification 

Particle swarm optimization (PSO) algorithm is simple and easy to implement. However, PSO can 
fall into the local optimum [14,15]. The original algorithm is modified and an adaptive escaping 
particle swarm optimization algorithm (AEPSO) is proposed. 

Supposing that X. = (x iV x i2 ,--,x id ) is the present position of ith particle; V t = (v n ,v. 2 ,---,v.^) is the 

present flying speed of ith particle; P i ={p i i->p i2 ->'"->Pid) * s the individual optimal fitness of ith 
particle; P g =(p g i,P g2 >'">P g d) * s g rou P optimal fitness; d is particle dimension, the evolution 
equations of original PSO algorithm are: 



where, k is iteration times, co is inertia weight; c x and c 2 are the random numbers between (0, 2) and 
called learning factors; r x and r 2 are the random numbers between (0, 1) [16,17]. 

When P g does not change over M generations, all particles are close to P g . The flying speed of 
particles is very little and subsistence density is large. An escape strategy should be adopted to update 
the velocity of particles and enlarge the searching space. The escape strategy changes the velocity of 
particles and makes variation in order to increase the diversity of particles and jump out of local 
optimization. When escaping, the particles are divided into two components, one component updates 



vf +1 = w\ + c x r x Of - x\ ) + c 2 r 2 (p k g - xf ) 



(14) 



(15) 



Sensors 2014, 14 



12693 



their velocity according to Equation (16) and the other component updates their velocity according to 
Equation (17). 

v k id +l =r,-v k id -{A-v k id y (16) 

V td =V^ max (17) 

where r 3 and r 4 are the random numbers between (0, 1). A is a constant that controls the particle 
velocity. v max is the maximal velocity value that particles have. 

The declining linearly inertia weight helps to swarm searching. The adjusting strategy is as follows: 

Hk) = w m ^- (w -~ w ^K k (18) 

^max 

where w(k) is the current inertia weight; w max and w min are the maximum and minimum inertia 
weights; q max is the maximum number of iterations; k is current iteration times. 

3.2.2. Optimizing Sensor Distribution Based on AEPSO Algorithm 

The adaptive escaping particle swarm optimization algorithm is adopted to obtain the optimum value of 
E and the sensor distribution. The detailed process is as follows and a flow chart is illustrated in Figure 3. 

(1) Producing a particle. A d-dimension particle is produced by selecting a distribution model 

randomly according to the number of sensors, which is expressed as 
a = {xr x ,'",xr m ,yr x ,--,yr m ), where xr j and yr. represent x and y coordinate of jth sensor, 

respectively, and m is the number of sensors; d = 2m . 

(2) Calculating matrix D and the rank of D. 

(3) If all column vectors of matrix D are not zero and matrix D is full rank, put this particle into 
the initial particle swarm; otherwise, return to Step (1). 

(4) Calculating the optimal group fitness when the number of particles in initial swarm is equal to 
the given number. 

(5) Updating the velocity and position of each particle according to Equations (14) and (15) and 
calculating fitness. 

(6) Updating the individual optimal fitness. If the current fitness is superior to experienced 
optimal fitness P. , the experienced optimal fitness P. is replaced by the current fitness. 

(7) Updating the group optimal fitness. If the current fitness is superior to the group optimal 
fitness, the group optimal fitness P g is replaced by current fitness. 

(8) Recording group optimal fitness P g in each iteration. If the value of P g does not change over 

M generations, return to Step (9) and adopt escaping strategy; otherwise, return to Step (5). 

(9) Dividing the particles into two components and updating particles velocity according to 
Equations (16) and (17). 

(10) If the ending condition is met, terminate the iteration or return to Step (5). 
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Figure 3. Process of optimizing sensor distribution. 
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4. Numerical Simulations 

4.1. Comparison of Cell Performance 

To verify the sensor distribution design, a test region is 16 m x 16 m, explosive is placed in the 
lower-left and sensors are distributed on the region boundary. Two cell patterns are used in dividing 
the test region. The first one is a uniform rectangle with 49 cells as shown in Figure 4a. The second 
one is sub-region and multi-scale cells as shown in Figure 4b. The velocity of shock wave decreases 
exponentially with the distance according to the prior information. In the near region to the explosive, 
the attenuation amplitude of velocity is larger and the cells are smaller and denser. In the far region the 
attenuation amplitude of velocity is smaller and the cells are bigger and sparser. The region is divided 
into seven sub-regions according to the proportional distance and different size cells in the symmetry 
region are placed in order to avoid a row vector to be linearity dependent on matrix D . The total cells 
number is 58. 




The influence of cells on matrix D is compared with the same sensor distribution. 

(1) To ensure each cell being passed through by at least one ray, the lowest number of sensor 
required by uniform rectangle cells is 9 as shown in Figure 4a. The multi-scale cells require at 
least five sensors as shown in Figure 4b. 

(2) Each index of matrix D in uniform rectangle cells and multi-scale cells is compared with 
13 sensors and the same distribution. Simulation results are given in Table 1. It is clear that the 
indexes in multi-scale cells are superior to those of the uniform rectangle cells. The matrix D in 
multi-scale cells is full rank mostly and the condition number is far smaller than that of 
uniform rectangle cells. The ray density and orthogonality in multi-scale cells are larger than 
that of uniform rectangle cells generally. 
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Table 1. Each index and simulation results. 



Cells Pattern 


Number 


Rank 


El 


E2 


O 


P 




1 


13 


53.08 


111.08 


0.220 


3.334 




2 


13 


55.38 


79.71 


0.225 


3.235 




3 


13 


55.57 


50.38 


0.229 


3.235 


Multi-Scale Cells 


4 


13 


56.04 


3285.90 


0.227 


3.353 




5 


13 


57.34 


1178.80 


0.185 


3.235 




0 




jy.jl 


oJ4.UU 


U.zUo 






7 

/ 




jy. 1 d 




U.Z lo 


'X zK1 




1 


12 


54.49 


1.76 x 10 


0.134 


2.653 




2 


12 


56.51 


6.89 x 10 15 


0.149 


2.755 


Uniform 
Rectangle Cells 


3 
4 
5 


12 
12 
12 


56.87 
57.58 
58.54 


1.02 x 10 16 
9.45 x 10 15 
8.86 x 10 15 


0.145 
0.151 
0.139 


2.694 
2.755 
2.857 




6 


11 


60.91 


1.24 x 10 16 


0.141 


2.857 




7 


11 


60.38 


2.14 x 10 16 


0.141 


2.898 



4.2. Optimizing Sensor Distribution Based on Multi-Scale Cells 

4.2.1. Parameters Setting and Simulation of AEPSO Algorithm 

To validate the AEPSO algorithm, two multi-modal functions are used to test the PSO and 
AEPSO algorithms. 
Rastrigrin function: 

/ 1 (x) = XU i 2 -10cos(2^.) + 10] (-10<x, <10) ( 19) 

Griewank function: 

& {x) = ^ s x i - n c ° s ^ 1 ^ + 1 ( - 5 °° * x > • * 5 °° ) (20) 

4UUU i= i i= i 

The two functions have many local minimum values, which has the global minimum value of 0 
when Xi = 0. The mean best fitness (MBF) and standard deviation (SD) are used to estimate the 
solutions. Parameters setting are given in Table 2. 



Table 2. Parameters setting. 



Algorithm 


CO 


c i 


c 2 


A 


M 


PSO 


0.7 


1.49 


1.49 






AEPSO 


[0.9, 0.4] 


1.49 


1.49 


1 


10 
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The simulation results are given in Table 3. From the table, it is concluded that AEPSO is superior 
to PSO algorithm in accuracy, convergence, and stability. 



Table 3. Simulation results. 



Function 


Dimension 


Algorithm 


MBF 


SD 


Rastrigrin Function 


10 


PSO 
AEPSO 


40.1 

1.8 x 10" 2 


27.4 
2.1 x 10 -2 


30 


PSO 


149.7 


45.3 




AEPSO 


8.5 x 10" 2 


1.4 x 10" 1 


Griewank Function 


10 


PSO 
AEPSO 


2.6 x 10" 1 
2.5 x 10 -4 


1.2 x 10" 1 
6.9 x 10" 4 


30 


PSO 
AEPSO 


6.8 
3.6 x 10" 2 


7.1 

4.5 x 10" 2 



4.2.2. Optimizing Sensor Distribution 

Each index based on multi-scale cells with symmetrical distribution and random distribution 
of sensors is shown in Figure 5a,b. According to the values of E x , E 2 , p and O , the weight 
coefficients are selected as: co x =0.8, co 2 =0.2, a> 3 =15, k x =0.1, k 2 =0.9. The adaptive escaping 
particle swarm optimization (PSO) algorithm is adopted to obtain the optimum value of E and 
sensor distribution as shown in Figure 5c. The values of E x , E 2 , p and O with different sensor 
distributions are in Table 4. It is clear that the indexes in optimum distribution are superior to that of 
the others: smallest E x andi? 2 , and the biggest E 3 . Figure 6 shows a comparison of the normalized 
eigenvalue spectra respectively in each case. The optimized model is again superior to that 
of the others. Numerical simulations are presented with the same initial model and inversion algorithm. 
From the results it can be concluded that the optimum distribution can result in the smallest 
inversion error. 



Table 4. Comparison of each index in different distributions. 



Sensor 


Symmetrical 


Random 


Optimum 


Distribution 


Distribution 


Distribution 


Distribution 


Rank 


11 


13 


13 


El 


56.68 


58.90 


50.84 


E2 


3.67 x 10 16 


887.70 


10.20 


E3 


0.48 


0.51 


0.57 


O 


0.178 


0.190 


0.254 


P 


3.157 


3.373 


3.373 


E 


7.35 x 10 15 


254.16 


69.24 


Related error (%) 


7.89 


7.41 


5.51 
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Figure 5. Different sensor distributions, respectively in each case, (a) Symmetrical sensor 
distribution; (b) Random sensor distribution; (c) Optimum sensor distribution. 
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5. Conclusions 

For travel time tomography in explosion, optimizing sensor distribution is very important for saving 
the costs and increasing the acquired information. This paper shows the inversion stability can be 
improved by designing optimal sensor distribution. This paper analyzed the effect of eigenvalue, rank, 
condition number and rays coverage on improving matrix D and proposes the evaluating function on 
optimizing sensor distribution. An adaptive escaping particle swarm optimization algorithm is used to 
obtain the optimum sensor distribution based on sub-region and multi-scale cells. Simulation shows 
that the optimization method is feasible and the optimal sensor design achieves inversion stability. 
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